#Clear R
rm(list=ls())

load("Kam_data.RData")


#function to install packages if they don't exist
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

# usage
packages <- c("melt", "ggplot","psych", "lavaan", "ks.test", "car", "lm.beta", "stargazer")
ipak(packages)

##############################
#
#Table A.1 Descriptive statistics
#
##############################

#DV
summary(data$DV_irradiation)
sd(data$DV_irradiation, na.rm=T)

#full 18 item
summary(data$NFC18)
sd(data$NFC18, na.rm=T)
alpha_NFC<-psych::alpha(data.frame(data$NfC_1, data$NfC_2, data$NfC_3_rec, data$NfC_4_rec, data$NfC_5_rec,data$NfC_6, data$NfC_7_rec, data$NfC_8_rec, data$NfC_9_rec,data$NfC_10, data$NfC_11, data$NfC_12_rec,data$NfC_13 , data$NfC_14 , data$NfC_15, data$NfC_16_rec,data$NfC_17_rec,data$NfC_18))
summary(alpha_NFC)

skewness(data$NFC18, na.rm=T)
kurtosis(data$NFC18, na.rm=T)
var(data$NFC18, na.rm=T)

#6-item measure
summary(data$NFC_Bullock)
sd(data$NFC_Bullock, na.rm=T)

alpha_NFC_Bullock<-psych::alpha(data.frame(data$NfC_1, data$NfC_2, data$NfC_3_rec, data$NfC_6, data$NfC_16_rec, data$NfC_8_rec))
summary(alpha_NFC_Bullock)
skewness(data$NFC_Bullock, na.rm=T)
kurtosis(data$NFC_Bullock, na.rm=T)
var(data$NFC_Bullock, na.rm=T)

#2-item measure
summary(data$NFC_ANES)
sd(data$NFC_ANES, na.rm=T)
alpha_NFC_ANES<-psych::alpha(data.frame(data$NfC_1, data$NfC_2))
summary(alpha_NFC_ANES)
skewness(data$NFC_ANES, na.rm=T)
kurtosis(data$NFC_ANES, na.rm=T)
var(data$NFC_ANES, na.rm=T)

#gender
table(data$female)

#Age
summary(data$age)
sd(data$age, na.rm=T)
#Race
table(data$race)
#calculate percentages
(709*100)/881
(93*100)/881
(5*100)/881
(39*100)/881
(1*100)/881
(34*100)/881

#Education
table(data$education)
(23*100)/881
(158*100)/881
(199*100)/881
(114*100)/881
(241*100)/881
(124*100)/881
(22*100)/881

#Income
summary(data$income)
sd(data$income, na.rm=T)

##############################
#
#Table A.2 Convergent validity of the NFC
#
##############################

cor(data$NFC18,data$NFC_Bullock,use="pairwise.complete.obs")
cor(data$NFC18,data$NFC_ANES,use="pairwise.complete.obs")
cor(data$NFC_Bullock,data$NFC_ANES,use="pairwise.complete.obs")

##############################
#
#Table A.3 Kolmornov Smirnov Tests
#
##############################

ks.test(data$NFC18,data$NFC_Bullock)
ks.test(data$NFC18,data$NFC_ANES)

##############################
#
#Figure A.1 Distributions of the Need for Cognition using different batteries
#
##############################

#NFC Histograms of the three variables
forggplot_NFC<- melt(data,measure.vars = c("NFC_ANES","NFC_Bullock","NFC18")) 
ggplot(forggplot_NFC,aes(x=value,colour=variable, fill=variable)) +  scale_colour_manual(values=c("red","chartreuse2", "blue")) + labs(title="", x="Need for Cognition", y="Frequency") +geom_histogram(, binwidth=.02, alpha=1, position=position_dodge()) + scale_fill_manual(values=c("red","chartreuse2", "blue"), labels=c("ANES 2-item", "Bullock 6-item","Cacioppo et al. 18-item"), guide=guide_legend(title = "Battery"))  + guides(colour=FALSE) + theme_bw()

##############################
#
#Table A.4 Confirmatory Factor Analysis NFC
#
##############################

#fit 
CFAModel1<-'NFC18 = ~ NA*NfC_1+ NfC_2+ NfC_3_rec+ NfC_4_rec+ NfC_5_rec+ NfC_6+ NfC_7_rec + NfC_8_rec + NfC_9_rec+ NfC_10+ NfC_11+NfC_12_rec+ NfC_13+ NfC_14+ NfC_15+ NfC_16_rec+ NfC_17_rec+ NfC_18
# fix variance of speed factor
NFC18 ~~ 1*NFC18'
#assess fit

fit<-cfa(CFAModel1, data=data)
summary(fit, fit.measures=TRUE, standardized=TRUE)
inspect(fit,what="std")

#*************************
#
#Table A.5 Study 1: Balance Checks different Need for Cognition Measures
#
#*************************

balance18<-(lm(data$NFC18~as.factor(data$Treatment),data=data))
balance6<-(lm(data$NFC_Bullock~as.factor(data$Treatment),data=data))
balanceANES2<-(lm(data$NFC_ANES~as.factor(data$Treatment),data=data))

#*************************
#
#Table A.6 Study 1: Balance Checks Background Characteristics
#
#*************************

#covariates
balanceSex<-glm(data$female~as.factor(data$Treatment), data=data, family=binomial(link="logit"))
balanceAge<-lm(data$age~as.factor(data$Treatment),data=data)
balanceEdu<-polr(as.factor(data$education)~as.factor(data$Treatment), data=data, Hess=TRUE)
balanceIncome<-lm(data$income~as.factor(data$Treatment),data=data)

#Extract results for models
stargazer(balanceSex, balanceAge, balanceEdu, balanceIncome, title="Balance Checks Demographics", align=TRUE, 
          dep.var.labels=c(""), 
          covariate.labels=c("Treatment"), 
          omit.stat=c("LL","ser","f", "adj.rsq"), 
          star.cutoffs=c(0.05), 
          notes = "OLS Regession models; *p<0.05", 
          notes.append = FALSE, 
          no.space=TRUE)

#*************************
#
#Table A.7 Study 1: Main Effects of Party Cues on Support for Food Irradiation
#
#*************************

MainEffect <-(lm(data$DV_irradiation~InParty +OutParty,data=data))
summary(MainEffect)

#*************************
#
#Table A.8 Study 1: Results Figure 1 from main text
#
#*************************
#model ANES
summary(NfCANES<-(lm(data$DV_irradiation~InParty*NFC+NoCues*NFC_ANES,data=data)))
#modle Bullock
summary(NfCBullock<-(lm(data$DV_irradiation~InParty*NFC_Bullock+NoCues*NFC_Bullock,data=data)))
#model full 18-item
summary(NFC18<-(lm(data$DV_irradiation~InParty*NFC18+NoCues*NFC18,data=data)))

#*************************
#
#Table A.9 Study 1: Standardized Results Figure 1 from main text
#
#*************************

ANES2 <-(lm(data$DV_irradiation~InParty*NFC_ANES +NoCues*NFC_ANES,data=data))
ANES2.beta<-lm.beta(ANES2)

Bullock6 <-(lm(data$DV_irradiation~InParty*NFC_Bullock +NoCues*NFC_Bullock,data=data))
Bullock6.beta<-lm.beta(Bullock6)

Full18 <-(lm(data$DV_irradiation~InParty*NFC18+NoCues*NFC18,data=data))
Full18.beta<-lm.beta(Full18)

#*************************
#
#Table A.10 Study 1: Structural equation model with the 18-item NfC measure
#
#*************************

data$OutIn[data$OutParty==1]=0
data$OutIn[data$InParty==1]=1

#create interaction variables
data$int1 <-(data$NfC_1*data$OutIn)
data$int2 <-(data$NfC_2*data$OutIn)
data$int3 <-(data$NfC_3_rec*data$OutIn)
data$int4 <-(data$NfC_4_rec*data$OutIn)
data$int5 <-(data$NfC_5_rec*data$OutIn)
data$int6 <-(data$NfC_6*data$OutIn)
data$int7 <-(data$NfC_7_rec*data$OutIn)
data$int8 <-(data$NfC_8_rec*data$OutIn)
data$int9 <-(data$NfC_9_rec*data$OutIn)
data$int10 <-(data$NfC_10*data$OutIn)
data$int11 <-(data$NfC_11*data$OutIn)
data$int12 <-(data$NfC_12_rec*data$OutIn)
data$int13 <-(data$NfC_13*data$OutIn)
data$int14 <-(data$NfC_14*data$OutIn)
data$int15 <-(data$NfC_15*data$OutIn)
data$int16 <-(data$NfC_16_rec*data$OutIn)
data$int17 <-(data$NfC_17_rec*data$OutIn)
data$int18 <-(data$NfC_18*data$OutIn)

#specify model
model.interaction18<-'
#latent variable
NFC18 = ~ NA*NfC_1+ NfC_2+ NfC_3_rec+ NfC_4_rec+ NfC_5_rec+ NfC_6+ NfC_7_rec+NfC_8_rec+NfC_9_rec+ NfC_10+ NfC_11+NfC_12_rec+ NfC_13+ NfC_14+ NfC_15+ NfC_16_rec+ NfC_17_rec+ NfC_18
Irradiation =~ SupportIrradiation_rec+ IrradiationGood_rec+ CostsFoodIrradiation_rec
Int =~ int1+int2+int3+int4+int5+int6+int7+int8+int9+int10+int11+int12+int13+int14+int15+int16+int17+int18
NFC18~~1*NFC18

#regression
Irradiation ~ NFC18+OutIn+Int'

#SEM model using Lavaan
fit_int18<-sem(model.interaction18 ,data=data, missing="listwise")
summary(fit_int18, fit.measures=TRUE, standardized=T)

################################################
#
#
#Appendix A.7 Study 1: Results using dependent variable employed by Kam (2005)
#
#
################################################

#recode variable to range from 0 - 1
data$SupportIrradiation_Kam<-((data$SupportIrradiation_rec-1)/4)

###########################################
#Make plot
############################################

#model ANES
summary(NfCANES<-(lm(data$SupportIrradiation_Kam~InParty*NFC_ANES +NoCues*NFC_ANES,data=data)))
#make a grpah
df1<-(interplot(m = NfCANES, var1 = "InParty", var2 = "NFC_ANES"))
df1<-(data.frame(df1$data))
df1$index<-"ANES 2-item"

#modle Bullock
summary(NfCBullock<-(lm(data$SupportIrradiation_Kam~InParty*NFC_Bullock +NoCues*NFC_Bullock,data=data)))
#make a graph
df2<-(interplot(m = NfCBullock, var1 = "InParty", var2 = "NFC_Bullock"))
df2<-(data.frame(df2$data))
df2$index<-"Bullock 6-item"

#model full 18-item
summary(NFC18<-(lm(data$SupportIrradiation_Kam~InParty*NFC18 +NoCues*NFC18,data=data)))
#make a graph
df3<-(interplot(m = NFC18, var1 = "InParty", var2 = "NFC18"))
df3<-(data.frame(df3$data))
df3$index<-"Cacioppo et al. 18-item"

toggplot<-(rbind(df1,df2,df3))
ggplot(toggplot,aes(x=fake,y=coef1))+geom_line()+geom_ribbon(aes(ymin=ub, ymax=lb), alpha=0.5)+facet_wrap(~index)+ geom_hline(yintercept=0, linetype = "dashed") + ylab('Marginal effect of In-Party Cue on Support for Ban on Food Irradiation') +xlab('Need for Cognition')+theme_bw()

################################################
#
#
#Appendix A.8 Study 1: Direct replication results Kam 2005
#
#
################################################

#Party Cue: 1 (strong or weak partisan, out-party endorses ban); 0.5 (leaning partisan, out- party endorses ban); 0 (no party cue); +0.5 (leaning partisan, in-party endorses ban); +1 (strong or weak partisan, in-party endorses ban).
data$PartyCue<-NA
data$PartyCue[data$pid1==1 & data$Treatment==2]<--1
data$PartyCue[data$pid1==2 & data$Treatment==3]<--1
data$PartyCue[data$pid4==1 & data$Treatment==2]<--0.5
data$PartyCue[data$pid4==2 & data$Treatment==3]<--0.5
data$PartyCue[data$Treatment==1 & data$pid1==1]<--0
data$PartyCue[data$Treatment==1 & data$pid1==2]<--0
data$PartyCue[data$Treatment==1 & data$pid4==1]<--0
data$PartyCue[data$Treatment==1 & data$pid4==2]<--0
data$PartyCue[data$Treatment==3 & data$pid4==1]<-0.5
data$PartyCue[data$Treatment==2 & data$pid4==2]<-0.5
data$PartyCue[data$Treatment==3 & data$pid1==1]<-1
data$PartyCue[data$Treatment==2 & data$pid1==2]<-1

#ANES 2-item
summary(NfCANES<-(lm(data$DV_irradiation~PartyCue*NFC_ANES,data=data)))

#Bullock 6-item
summary(NfCBullock<-(lm(data$DV_irradiation~PartyCue*NFC_Bullock,data=data)))

#NFC 18-item
summary(NFC18<-(lm(data$DV_irradiation~PartyCue*NFC18,data=data)))

stargazer(NfCANES, NfCBullock, NFC18, title="Moderation of Party Cues Following Kam (2005)", align=TRUE, 
          dep.var.labels=c("Support for Ban of Food Irradiation"), 
          covariate.labels=c("In-Party Cue", "Need for Cognition", "No-Party Cue", "In-Paty Cue X Need for Congition", "No Cue X Need for Cognition"), 
          omit.stat=c("LL","ser","f", "adj.rsq"), 
          star.cutoffs=c(0.05), 
          notes = "OLS Regession models; *p<0.05", 
          notes.append = FALSE, 
          no.space=TRUE)

